{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Glucose Minimal Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous chapter we implemented the glucose minimal model using given parameters, but I didn't say where those parameters came from.\n",
    "This notebook solves the mystery.\n",
    "\n",
    "We'll use a SciPy function called `leastsq`, which stands for \"least squares\"; that is, it finds the parameters that minimize the sum of squared differences between the results of the model and the data.\n",
    "\n",
    "You can think of `leastsq` as optional material.  We won't use it in the book itself, but it appears in a few of the case studies.\n",
    "\n",
    "You can read more about `leastsq` in [the SciPy documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html).  It uses the Levenberg-Marquart algorithm, which you can read about [on Wikipedia](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells download and read the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ModSim/raw/main/data/' +\n",
    "         'glucose_insulin.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('glucose_insulin.csv', index_col='time');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use `make_system` and `slope_func` as defined in Chapter 18."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
    "         'chap18.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from chap18 import make_system\n",
    "from chap18 import slope_func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing errors\n",
    "\n",
    "In this context, the \"errors\" are the differences between the results from the model and the data.\n",
    "\n",
    "To compute the errors, I'll start again with the parameters we used in Chapter 18."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "G0 = 270\n",
    "k1 = 0.02\n",
    "k2 = 0.02\n",
    "k3 = 1.5e-05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = G0, k1, k2, k3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system` takes the parameters and actual data and returns a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run the ODE solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_solve_ivp(system, slope_func, \n",
    "                                 t_eval=data.index)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we specify `t_eval=data.index`, the results are evaluated at the some time stamps as the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the results like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.glucose.plot(style='o', alpha=0.5, label='data')\n",
    "results.G.plot(style='-', color='C0', label='model')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration (mg/dL)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the first three time steps, the model does not fit the data. That's because it takes some time for the injected glucose to disperse.\n",
    "\n",
    "We can compute the errors like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "errors = results.G - data.glucose\n",
    "errors.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first few errors are substantially larger than the rest.\n",
    "\n",
    "In the next section, we'll use `leastsq` to search for the parameters that minimize the sum of the squared errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization\n",
    "\n",
    "To use `leastsq`, we need an \"error function\" that takes a sequence of parameters and returns an array of errors.\n",
    "\n",
    "Here's a function that does what we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_func(params, data):\n",
    "    \"\"\"Computes an array of errors to be minimized.\n",
    "    \n",
    "    params: sequence of parameters\n",
    "    data: DataFrame of values to be matched\n",
    "    \n",
    "    returns: array of errors\n",
    "    \"\"\"\n",
    "    print(params)\n",
    "    \n",
    "    # make a System with the given parameters\n",
    "    system = make_system(params, data)\n",
    "    \n",
    "    # solve the ODE\n",
    "    results, details = run_solve_ivp(system, slope_func, \n",
    "                                     t_eval=data.index)\n",
    "    \n",
    "    # compute the difference between the model\n",
    "    # results and actual data\n",
    "    errors = results.G - data.glucose\n",
    "    return errors.iloc[3:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`error_func` uses the given parameters to make a `System` object, runs the simulation, and returns the errors.\n",
    "\n",
    "But notice that it does not return all of the errors; rather, it uses `iloc` to select only the elements with index 3 or more.\n",
    "In other words, it omits the elements with index 0, 1, and 2.\n",
    "Note: You can read more about this use of `iloc` [in the Pandas documentation](https://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-integer).\n",
    "\n",
    "Since we don't expect the model to fit the data in this regime, we'll leave it out.\n",
    "\n",
    "We can call `error_func` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_func(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to call `leastsq`.  As arguments, we pass `error_func`, the parameters where we want to start the search, and the data, which will be passed as an argument to `error_func`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params, fit_details = leastsq(error_func, params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each time `error_func` is called, it prints the parameters, so we can get a sense of how `leastsq` works.\n",
    "\n",
    "`leastsq` has two return values.\n",
    "The first is an array with the best parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second is an object with information about the results, including a success flag and a message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_details.success"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_details.mesg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have `best_params`, we can use it to make a `System` object and run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "system2 = make_system(best_params, data)\n",
    "results2, details = run_solve_ivp(system2, slope_func, \n",
    "                                  t_eval=data.index)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results, along with the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.glucose.plot(style='o', alpha=0.5, label='data')\n",
    "results.G.plot(style='-', color='C0', label='model')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration (mg/dL)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the errors like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "errors2 = results2.G - data.glucose\n",
    "errors2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the sum of the squared errors like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import sum\n",
    "\n",
    "sum(errors2.iloc[3:]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If things have gone according to plan, the sum of squared errors should be smaller now, compared to the parameters we started with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(errors.iloc[3:]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting parameters\n",
    "\n",
    "So we found the parameters that best match the data.  You might wonder why this is useful.\n",
    "\n",
    "The parameters themselves don't mean very much, but we can use them to compute two quantities we care about:\n",
    "\n",
    "-   \"Glucose effectiveness\", $E$, which is the tendency of elevated\n",
    "    glucose to cause depletion of glucose.\n",
    "\n",
    "-   \"Insulin sensitivity\", $S$, which is the ability of elevated blood\n",
    "    insulin to enhance glucose effectiveness.\n",
    "\n",
    "Glucose effectiveness is defined as the change in $dG/dt$ as we vary\n",
    "$G$: \n",
    "\n",
    "$$E \\equiv - \\frac{\\delta \\dot{G}}{\\delta G}$$ \n",
    "\n",
    "where $\\dot{G}$ is shorthand for $dG/dt$. Taking the derivative of $dG/dt$ with respect to $G$, we get \n",
    "\n",
    "$$E = k_1 + X$$ \n",
    "\n",
    "The **glucose effectiveness index**, $S_G$, is the value of $E$ when blood insulin is near its basal level, $I_b$.\n",
    "In that case, $X$ approaches 0 and $E$ approaches $k_1$. So we can use\n",
    "the best-fit value of $k_1$ as an estimate of $S_G$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Insulin sensitivity is defined as the change in $E$ as we vary $I$:\n",
    "\n",
    "$$S \\equiv - \\frac{\\delta E}{\\delta I}$$ \n",
    "\n",
    "The **insulin sensitivity index**, $S_I$, is the value of $S$ when $E$ and $I$ are at steady state: \n",
    "\n",
    "$$S_I \\equiv \\frac{\\delta E_{SS}}{\\delta I_{SS}}$$ \n",
    "\n",
    "$E$ and $I$ are at steady state when $dG/dt$ and $dX/dt$ are 0, but we don't actually have to solve those equations to find $S_I$. \n",
    "\n",
    "If we set $dX/dt = 0$ and solve for $X$, we find the relation:\n",
    "\n",
    "$$X_{SS} = \\frac{k_3}{k_2} I_{SS}$$ \n",
    "\n",
    "And since $E = k_1 + X$, we have:\n",
    "\n",
    "$$S_I = \\frac{\\delta E_{SS}}{\\delta I_{SS}} = \\frac{\\delta X_{SS}}{\\delta I_{SS}}$$\n",
    "\n",
    "Taking the derivative of $X_{SS}$ with respect to $I_{SS}$, we have:\n",
    "\n",
    "$$S_I = k_3 / k_2$$ \n",
    "\n",
    "So if we find parameters that make the model fit the data, we can use $k_3 / k_2$ as an estimate of $S_I$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the parameters we found, here are the estimated values of $S_G$ and $S_I$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "G0, k1, k2, k3 = best_params\n",
    "indices = SimpleNamespace(S_G=k1, S_I=k3/k2)\n",
    "show(indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to [Boston et al](https://www.researchgate.net/publication/8931437_MINMOD_Millennium_A_Computer_Program_to_Calculate_Glucose_Effectiveness_and_Insulin_Sensitivity_From_the_Frequently_Sampled_Intravenous_Glucose_Tolerance_Test), normal ranges for these values are..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "S_G_interval = 1.2e-3, 4.5e-2\n",
    "S_G_interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "S_I_interval = 5.0e-5, 2.2e-3\n",
    "S_I_interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated quantities are within the normal intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:** How sensitive are the results to the starting guess for the parameters?  If you try different values for the starting guess, do we get the same values for the best parameters?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
